Load libraries

library(tidyverse)
library(scran)
library(scater)
library(gprofiler2)
library(gridExtra)
library(ComplexHeatmap)
library(circlize)
library(edgeR)
library(readxl)
library(ggrepel)
library(tftargets)
library(network)
library(ggnet)
library(parallel)
library(BiocParallel)
library(ggsci)
library(airr)
library(numbat)
library(scales)

ncores <- 10
mcparam <- MulticoreParam(workers=ncores)
register(mcparam)

Define variables

opt <- list()
opt$pathinfo <- "misc/pathways_2.xlsx"
opt$scebulk <- "data/submission/RNA_complete.RDS"
opt$plot <- "plots/Fig3/"


## Set theme
lgd <-  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5), 
        text = element_text(size = 15),
        panel.background = element_rect(fill = "transparent",colour = NA),
        plot.background = element_rect(fill = "transparent",colour = NA), 
        legend.background = element_rect(fill='transparent',colour = NA), 
        strip.background = element_blank(), 
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

Load data

## SCE
sce <- readRDS(opt$scebulk)
sce
## class: SingleCellExperiment 
## dim: 36601 117 
## metadata(0):
## assays(2): counts logcounts
## rownames(36601): ENSG00000000003 ENSG00000000005 ... ENSG00000288459
##   ENSG00000288460
## rowData names(3): ID Symbol Type
## colnames(117): H358_BafilomycinA1 H358_Birinapant ... H524_Selinexor
##   H524_Venetoclax
## colData names(5): uniqueBarcode patientID Diagnosis_simple Treatment
##   label
## reducedDimNames(2): PCA UMAP
## mainExpName: NULL
## altExpNames(0):
## Pathway information
path.info <- read_excel(opt$pathinfo)

Analysis

Perform differential gene expression testing

testTreat <- unique(sce$Treatment)

de.res <- mclapply("DMSO", mc.cores = ncores, 
                 function(z) {
                   
  message(paste0("Fitting model for ", z))
                   
  compDrug <- testTreat[testTreat != z] # instead of subtypes we use other treatments here
                   
  mclapply(compDrug, mc.cores = ncores, 
         function(x) {
    message(paste0("Drug ", x))
    summed <- sce#[hvgenes, ]
    colData(summed)$Treatment <- factor(colData(summed)$Treatment, levels = c(z, compDrug))
    y <- DGEList(counts(summed), samples = colData(summed))
    design <- model.matrix(~0 + factor(patientID) + Treatment, y$samples) # set coefficient (name or matrix column) - fig. 7 f1000
    keep <- filterByExpr(y, design = design)
    y <- y[keep,]
    summary(keep)
    message(paste0("Keeping ", sum(keep), " genes"))
    y <- calcNormFactors(y)
    y <- estimateDisp(y, design) # gene-wise dispersion (variance on top of poisson distribution)
    plotBCV(y)
    fit <- glmQLFit(y, design, robust = TRUE) # QL way of fitting the model
    plotQLDisp(fit)
    res <- glmQLFTest(fit, coef = paste0("Treatment", x))
    summary(decideTests(res))
    resTab <- res$table
    resTab$gene <- rownames(resTab)
    resTab <- resTab %>% mutate(p.adj = p.adjust(PValue, method = "BH")) #%>%
      #mutate(p.adj = ifelse(p.adj < 1e-12, 1e-12, p.adj))
    resTab$Treatment <- x
    resTab$Control <- z
    resTab %>% arrange(desc(logFC))
  }) %>% bind_rows()
}) %>% bind_rows()
## Fitting model for DMSO
de.res <- de.res %>%
  left_join(select(data.frame(rowData(sce)), ID, Symbol), c("gene" = "ID"))

de.res %>% filter(Treatment == "Selinexor", Symbol %in% c("BCL2", "TP53"))
##        logFC   logCPM           F       PValue            gene        p.adj
## 1 0.99175500 6.148050 148.6069142 8.633760e-22 ENSG00000141510 6.892299e-19
## 2 0.07380089 6.638112   0.2161808 6.429888e-01 ENSG00000171791 8.493301e-01
##   Treatment Control Symbol
## 1 Selinexor    DMSO   TP53
## 2 Selinexor    DMSO   BCL2
#write.csv(de.res, opt$deres)

Show the number of DEG

sigTab <- de.res %>%
  group_by(Treatment) %>%
  filter(p.adj < 0.1) %>%
  filter(logFC > 0.25 | logFC < -0.25) %>%
  mutate(dir = ifelse(logFC > 0, "Up", "Down")) %>%
  dplyr::count(dir) %>%
  dplyr::rename(Direction = dir) %>%
  ungroup()

drugOrder <- sigTab %>%
  group_by(Treatment) %>%
  summarise(allDEG = sum(n)) %>%
  arrange(desc(allDEG)) %>%
  pull(Treatment)

sigTab <- sigTab %>%
  mutate(Treatment = factor(Treatment, levels = drugOrder), 
         n = ifelse(Direction == "Down", n * -1, n))

p <- sigTab %>%
  ggplot(aes(x = Treatment, y = n, fill = Direction)) +
  geom_bar(stat = "identity", colour = "black") +
  geom_text(aes(label = n), data = sigTab[sigTab$Direction == "Up", ], size = 5, nudge_y = 200) +
  geom_text(aes(label = n), data = sigTab[sigTab$Direction == "Down", ], size = 5, nudge_y = -200) +
  xlab("") + ylab("Number of Genes") +
  ggtitle("Significantly Differentially Expressed Genes") +
  lgd + scale_fill_manual(values = c("Up" = "#F44336", "Down" = "#1976D2")) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1), 
        plot.title = element_text(size = 21, vjust = 1), 
        axis.text = element_text(size = 17.5), 
        text = element_text(size = 22.5), 
        panel.border = element_rect(linewidth = 1))
p

ggsave(plot = p, filename = paste0(opt$plot, "n_deg.png"), 
       height = 17.5, width = 22, units = "cm")

Construct a network of differentially expressed genes

# Assign nodes (treatments) and edges (5% FDR DEGs)
plotMat <- de.res %>% filter(p.adj < 0.1) %>%
  filter(logFC > 0.25 | logFC < -0.25)
allTreat <- unique(plotMat$Treatment)
target <- data.frame(Treatment = allTreat,
                     target = c("Autophagyi", "IAPi", "HDACi", "mTORi", "ITKi", "TLR8a", "MDM2i", 
                                "IMID", "JAKi", "XPOi", "BCL2i"))

opt$sub <- FALSE

if(opt$sub == TRUE) {
  plotMat <- de.res %>% filter(p.adj < 0.1, !Treatment %in% c("Pomalidomide", "Dacinostat")
                             ) %>%
  filter(logFC > 0.25 | logFC < -0.25) #%>%
    #filter(logFC < -0.25)
} else {
  plotMat <- de.res %>% filter(p.adj < 0.1) %>%
    filter(logFC > 0.25 | logFC < -0.25) #%>%
    #filter(logFC > 0.25)
}


allTreat <- unique(plotMat$Treatment)

nodes <- data.frame(Treatment = allTreat)
edges <- plotMat %>% select(Treatment, Symbol) %>%
  mutate(Treatment = as.character(Treatment))

# Create network object
set.seed(100)
net <- network(edges, directed = TRUE)
netTab <- ggnet2(net)$data |> 
  mutate(Type = ifelse(label %in% allTreat, "Treatment", "Symbol"),
         Treatment = ifelse(Type == "Treatment", label, NA)) |> 
  select(-alpha,-color,-shape,-size) |> 
  mutate(x = as.numeric(x), y = as.numeric(y))
  
# Define edges based on direction and padj of DEGs
from <- edges %>%
  left_join(netTab, by = c("Treatment" = "label")) %>% dplyr::select(Treatment, x, y) %>%
  dplyr::rename(x_start = x, y_start = y) |> distinct()

to <- edges %>%
  left_join(netTab, by = c("Symbol" = "label")) %>% dplyr::select(Symbol,x,y) %>%
  dplyr::rename(x_end = x, y_end = y) |> distinct()

edges <- edges %>% left_join(from) %>% left_join(to) %>%
  left_join(plotMat %>% select(Symbol,Treatment,logFC,p.adj)) %>%
  distinct() %>%
  mutate(Direction = ifelse(logFC < 0,"Down","Up"),
         alpha = -log10(p.adj)
         #alpha = rescale(-log10(p.adj#squish(p.adj, c(1.812407e-209,1))))
         )

# Create auxiliary table for annotation

n_DEG <- plotMat %>% #filter(!Treatment %in% c("Pomalidomide", "BafilomycinA1", "Dacinostat")) %>%
  group_by(Treatment) %>%
  filter(p.adj < 0.1) %>%
  filter(logFC > 0.25 | logFC < -0.25) %>% 
  dplyr::count() %>% ungroup()

netTabAnno <- netTab |> 
  left_join(target) |> 
  left_join(n_DEG) |> 
  mutate(Treatment = factor(Treatment, levels = allTreat))

p <- ggplot()  +
  geom_curve(data = edges,
               aes(x = x_start, xend = x_end, y = y_start, yend = y_end, col = Direction, alpha = p.adj),
               linewidth = 0.1, curvature = 0.05) +
  geom_point(data = netTabAnno[netTabAnno$Type == "Symbol",],
             aes(x = x, y = y),
             size = 0.5, shape = 21, color = "grey60", fill = "grey95", stroke = 0.75, alpha = 0.6) +
  geom_point(data = netTabAnno[netTabAnno$Type == "Treatment",],
             aes(x = x, y = y, fill = Treatment, size = n), shape = 21, stroke = 0.75) +
  # scale_fill_manual(values = col_TX) +
  scale_size_continuous(name = "Number of DEGs",
                        breaks = seq(0, max(n_DEG$n), by = 1000), limits = c(0, max(n_DEG$n)), range = c(0.5,7.5)) +
  scale_color_manual(values = c("Up" = "#F44336", "Down" = "#1976D2")) +
  scale_alpha_continuous() +
  geom_label_repel(data = netTabAnno |> filter(Type == "Treatment"),
                  aes(x = x, y = y, label = Treatment, fill = Treatment, point.size = rescale(n, to = c(0.5,7.5))),
                  max.overlaps = Inf, vjust = "bottom", nudge_y = -0.001, fontface = "bold", alpha = 0.7, #direction = "y", 
                  show.legend = FALSE, size = 5)  +
  guides(color = guide_legend(ncol = 2, override.aes = list(linewidth=2.5)),
         size = guide_legend(ncol = 2),
         alpha = guide_legend(override.aes = list(linewidth=2.5)),
         fill = guide_legend(override.aes = list(size=3.5))) +
  xlab("") + ylab("") + ggtitle("Network of Differentially Expressed Genes") +
  lgd +
  theme(axis.ticks = element_blank(), 
        axis.text = element_blank(), 
        panel.border = element_rect(linewidth = 1, colour = NA, fill = "transparent")) +
  scale_fill_manual(values = c("Venetoclax" = "grey70", "Birinapant" = "#AD1457", "Nutlin-3a" = "#FFA000", 
                                "Ibrutinib" = "#283593", "Everolimus" = "#43A047", "Selinexor" = "#FFD45F",
                               "BafilomycinA1" = "#F8BBD0", "Motolimod" = "#AB47BC", "Dacinostat" = "#F44336", 
                                "Dacinostat" = "#64B5F6", "Pomalidomide" = "#80DEEA", "Ruxolitinib" = "#4DB6AC"))

p

if(opt$sub == TRUE) {
 ggsave(plot = p, filename = paste0(opt$plot, "deg_subnetwork.png"),
       height = 20, width = 25, units = "cm")
} else {
  ggsave(plot = p, filename = paste0(opt$plot, "deg_network.png"),
       height = 20, width = 25, units = "cm")
}

Volcano plots

pList <- lapply(c("Birinapant", "Selinexor", "Nutlin-3a"), function(x) {
  
  if(x %in% c("Motolimod", "Everolimus", "Ibrutinib", "BaflimoycinA1")) {
    geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()
  } else if(x %in% c("Birinapant")) {
    geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF")) %>% select(gene) %>% unlist() %>% as.vector()
  } else if(x %in% c("Nutlin-3a", "Selinexor")) {
    geneSelec <- path.info %>% filter(pathway %in% c("TP53up")) %>% select(gene) %>% unlist() %>% as.vector()
    geneSelec <- c("TP53", geneSelec)

  } else if(x %in% c("Everolimus")) {
    geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()
  }
  
  volTab <- de.res %>% filter(Treatment == x)
  #volTab <- de.res.2020 %>% filter(Treatment == x, celltype == "T-PLL") %>% dplyr::rename(Symbol = gene)
  
  ## Define limits  
  minLim <- -1 * max(abs(de.res$logFC))
  maxLim <- max(abs(de.res$logFC))
  thresh <- 0.25

  volTab$diffCol <- NA
  volTab[volTab$logFC > 0 & volTab$p.adj < 0.1 & volTab$logFC > thresh, "diffCol"] <- "sigUp"
  volTab[volTab$logFC < 0 & volTab$p.adj < 0.1 & volTab$logFC < -thresh, "diffCol"] <- "sigDown"
  volTab[volTab$p.adj >= 0.1, "diffCol"] <- "ns"
  volTab[volTab$logFC < thresh & volTab$logFC > -(thresh), "diffCol"] <- "ns"
  
  p <- volTab %>%
    mutate(logFC = ifelse(logFC > maxLim, maxLim, logFC)) %>%
    mutate(logFC = ifelse(logFC < minLim, minLim, logFC)) %>%
    #filter(gene %in% geneSelec) %>%
    filter(Treatment == x) %>%
    ggplot(aes(x = logFC, y = -log10(p.adj), fill = diffCol)) +
    geom_point(shape = 21, size = 2) + xlab("Log2FC") + ylab("-Log10(p.adj)") + ggtitle(paste0(x)) +
    #xlim(minLim, maxLim) +
    geom_hline(yintercept = -log10(0.1), linetype = "dashed") +
    geom_vline(xintercept = -thresh, linetype = "dashed") + geom_vline(xintercept = thresh, linetype = "dashed") +
    geom_label_repel(data = filter(volTab, Symbol %in% c(geneSelec), Treatment %in% x, p.adj < 0.1, 
                                 logFC > thresh | logFC < -(thresh)), aes(label = Symbol), fill = "white", 
                    nudge_y = 0.25, size = 4, max.overlaps = getOption("ggrepel.max.overlaps", default = 30), 
                     segment.linetype = 2, force = 20, alpha = 0.8) +
    scale_fill_manual(values = c("sigUp" = "#F44336", "sigDown" = "#1976D2", "ns" = "grey")) +
    theme_classic() +
    theme(plot.title = element_text(hjust = 0.5), 
          legend.position = "none", 
          text = element_text(size = 17.5), 
          axis.text = element_text(size = 15),
          panel.background = element_rect(fill = "transparent",colour = NA),
          plot.background = element_rect(fill = "transparent",colour = NA),
          legend.background = element_rect(fill='transparent'),
          strip.background = element_blank())
  
   ggsave(plot = p, filename = paste0(opt$plot, "vol_", x, ".png"),
          height = 15, width = 15, unit = "cm")
  p

})
## Warning: ggrepel: 12 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
## Warning: ggrepel: 9 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps
pList 
## [[1]]

## 
## [[2]]
## Warning: ggrepel: 14 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

## 
## [[3]]
## Warning: ggrepel: 10 unlabeled data points (too many overlaps). Consider
## increasing max.overlaps

geneSelec <- path.info %>% filter(pathway %in% c("NFKB", "NFKB_Chemo", "NFKB_GF", "NFKB_TF", "NFKB_Apo", "TNF", "TLR")) %>% select(gene) %>% unlist() %>% as.vector()

geneSub <- lapply(unique(de.res$Treatment), function(x) {
  de.res %>% 
    filter(Treatment == x) %>%
    filter(p.adj < 0.1) %>%
    filter(logFC > 0.25 | logFC < -0.25) %>%
    pull(Symbol)
  
}) %>% unlist() %>% unique()

geneSub <- intersect(geneSelec, geneSub)

plotMat.wide <- de.res %>% 
  mutate(logFC = ifelse(p.adj >= 0.1, 0, logFC)) %>%
  dplyr::select(Symbol, Treatment, logFC) %>% 
  filter(Symbol %in% geneSub) %>%
  dplyr::rename(log2FC = logFC) %>%
  pivot_wider(names_from = "Treatment", values_from = "log2FC") %>%
  column_to_rownames("Symbol")

geneOrder <- rownames(plotMat.wide)[hclust(dist(plotMat.wide), method = "ward.D2")$order]

Dot plot

p <- de.res %>% filter(#Treatment %in% c("Birinapant", "Motolimod"),
  Symbol %in% geneSub) %>%
  filter(p.adj < 0.1) %>%
  mutate(pSign = -log10(p.adj) * sign(logFC),
         p.adj = ifelse(-log10(p.adj) > 5, 5, -log10(p.adj)), 
         Symbol = factor(Symbol, levels = geneOrder)) %>%
  dplyr::rename(log2FC = logFC) %>%
  ggplot(aes(x = Symbol, y = Treatment)) +
  geom_point(aes(size = p.adj, fill = log2FC), shape = 21) +
  xlab("") + ylab("") + ggtitle("Effect on TNF-NFkB and TLR Target Genes") +
  lgd +
  scale_fill_gradient2(high = "#F44336", low = "#1976D2", name = "log2FC",
                       limits=c(-2, 2), oob=squish) +
  scale_size_continuous(name = "-Log10(p.adj)",
                        breaks = seq(0, 5, by = 2.5),
                        limits = c(0, 5), #range = c(0.5,7.5)
                        ) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1),
        text = element_text(size = 12.5),
        panel.border = element_rect(linewidth = 1))
p

ggsave(plot = p, filename = paste0(opt$plot, "biri_moto_nfkb.png"),
       height = 9, width = 50, units = "cm")

GSEA

go.df.merge <- mclapply(c("Birinapant", "Motolimod", "Selinexor", "BafilomycinA1", "Dacinostat", "Pomalidomide", "Ibrutinib", 
                        "Everolimus", "Venetoclax", "Ruxolitinib", "Nutlin-3a"), mc.cores = ncores, 
                      function(compTreat) {
    message(paste("Running gprofiler for", compTreat, sep = " "))
    sigUp <- de.res %>% filter(logFC > 0.25, Treatment == compTreat, 
                               #!gene %in% c(drugOther), 
                               p.adj <= 0.1) %>% arrange(desc(logFC))
    #sigUp <- de.res %>% filter(celltype == compType, logFC > 0.25, Treatment == compTreat, p.adj <= 0.05) %>% arrange(PValue)
    
    sigUp
    if(nrow(sigUp) > 1) {
      gostresUp <- gost(query = sigUp$gene,
                        organism = "hsapiens", ordered_query = TRUE,
                        multi_query = FALSE, significant = TRUE, exclude_iea = FALSE,
                        measure_underrepresentation = FALSE, evcodes = FALSE,
                        user_threshold = 0.1, correction_method = "fdr",
                        domain_scope = "annotated", custom_bg = NULL,
                        numeric_ns = "", sources = NULL, as_short_link = FALSE)
      gostresUp$result

      go.df.up <- gostresUp$result %>% arrange(p_value)
      go.df.up$rank <- 1:nrow(go.df.up)
      go.df.up$type <- "Up"
    } else {
      go.df.up <- data.frame(p_value = NA, rank = NA, type = "Up", source = NA, term_name = NA, 
                               intersection_size = 10)
    }

    sigDown <- de.res %>% filter(logFC < -(0.25), Treatment == compTreat,
                                 #!gene %in% c(drugOther),
                                  p.adj <= 0.1) %>% arrange(logFC)
    #sigDown <- de.res %>% filter(celltype == compType, logFC < -(0.25), Treatment == compTreat, p.adj <= 0.05) %>% arrange(PValue)
    if(nrow(sigDown) > 1) {
    gostresDown <- gost(query = sigDown$gene,
                        organism = "hsapiens", ordered_query = TRUE,
                        multi_query = FALSE, significant = TRUE, exclude_iea = FALSE,
                        measure_underrepresentation = FALSE, evcodes = FALSE,
                        user_threshold = 0.1, correction_method = "fdr",
                        domain_scope = "annotated", custom_bg = NULL,
                        numeric_ns = "", sources = NULL, as_short_link = FALSE)

    gostresDown$result

    go.df.down <- gostresDown$result %>% arrange(p_value)
    go.df.down$rank <- 1:nrow(go.df.down) * -1
    go.df.down$type <- "Down"
    } else {
      go.df.down <- data.frame(p_value = NA, rank = NA, type = "Down", source = NA, term_name = NA, 
                               intersection_size = 10)
    }

    go.df.merge <- rbind(go.df.up, go.df.down)
    go.df.merge$Treatment <- compTreat

    go.df.merge <- go.df.merge %>% filter(!intersection_size < 10) %>% 
      filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
      mutate(p.adj = p.adjust(p_value, method = "BH"))
    
    if(nrow(go.df.merge) < 1) {
      go.df.merge <- data.frame(p_value = NA, rank = NA, type = "Down", source = NA, term_name = NA, 
                               intersection_size = 0)
    } else {
      go.df.merge
    }
}) %>% bind_rows() 

Show number of enriched gene-sets

#### Show number of enriched gene sets
sigTab <- go.df.merge %>%
  dplyr::count(Treatment, type) %>%
  dplyr::rename(Direction = type) %>%
  mutate(n = ifelse(Treatment %in% c("Ruxolitinib", "Venetoclax"), 0, n)) %>%
  mutate(Treatment = factor(Treatment, levels = drugOrder), 
         n = ifelse(Direction == "Down", n * -1, n)) %>%
  filter(!is.na(Treatment))

naTab <- data.frame(Treatment = c("Venetoclax", "Venetoclax", "Selinexor"), 
                    Direction = c("Down", "Up", "Down"), 
                    n = c(0, 0, 0))
sigTab <- rbind(sigTab, naTab)

p <- sigTab %>%
  ggplot(aes(x = Treatment, y = n, fill = Direction)) +
  geom_bar(stat = "identity", colour = "black") +
  geom_text(aes(label = n), data = sigTab[sigTab$Direction == "Up", ], size = 5, nudge_y = 200) +
  geom_text(aes(label = n), data = sigTab[sigTab$Direction == "Down", ], size = 5, nudge_y = -200) +
  xlab("") + ylab("Number of Gene Sets") +
  ggtitle("Significantly Enriched Gene Sets") +
  lgd + scale_fill_manual(values = c("Up" = "#F44336", "Down" = "#1976D2")) +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1), 
        plot.title = element_text(size = 21, vjust = 1), 
        axis.text = element_text(size = 17.5), 
        text = element_text(size = 22.5), 
        panel.border = element_rect(linewidth = 1))
p

ggsave(plot = p, filename = paste0(opt$plot, "n_gsea.png"), 
       height = 17.5, width = 22, units = "cm")

Make dot plot

## Exclude all terms that are up and downregulated by a single-drug
dblTerm <- go.df.merge %>% 
  filter(!intersection_size < 10) %>%
  filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
  mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
  #filter(term_name %in% unique(termTab$term_name)) %>%
  dplyr::count(term_name, Treatment) %>%
  filter(n > 1) %>%
  pull(term_name) %>% unique()

## Keep the top20 of each compound
nTop <- 5
termTab <- lapply(c("Selinexor", "Nutlin-3a", "Venetoclax", "Birinapant", "Ibrutinib", "Motolimod", 
                    "Dacinostat", "Pomalidomide", "Everolimus"), function(y) {
  termTab <- lapply(c("Up", "Down"), function(z) {
    message(y)
    message(z)
      
    gseaTab <- go.df.merge %>% filter(type == z, Treatment == y) %>% 
      filter(!intersection_size < 10) %>%
      filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
      mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
      filter(!term_name %in% dblTerm)
      
    gseaTab <- gseaTab[1:nTop, ] # get only top enrichments
      
    #termOrder <- gseaTab %>% arrange(desc(p_value)) %>% pull(term_name)
      
    gseaTab <- gseaTab %>%
      mutate(p.adj = p.adjust(p_value, method = "BH")) 
  }) %>% bind_rows()
}) %>% bind_rows()

gseaTab <- go.df.merge %>% 
  filter(!intersection_size < 10) %>%
  filter(!source %in% c("CORUM", "TF", "MIRNA")) %>%
  mutate(term_name = paste0(source, "_", term_name)) %>% unique() %>%
  filter(term_name %in% unique(termTab$term_name)) %>%
  filter(!term_name %in% dblTerm)

plotMat <- gseaTab %>%

    mutate(p.adj = -log10(p.adj), p.adj = ifelse(type == "Down", p.adj * -1, p.adj)) %>%
    select(term_name, Treatment, p.adj) %>% unique() %>%
    pivot_wider(names_from = "Treatment", values_from = "p.adj") %>%
    column_to_rownames("term_name")

plotMat[is.na(plotMat)] <- 0
  
geneOrder <- rownames(plotMat)[hclust(dist(plotMat), method = "ward.D2")$order]

p <- gseaTab %>%
  mutate(term_name = factor(term_name, levels = geneOrder), 
         intersection_size = ifelse(intersection_size > 200, 200, intersection_size)) %>%
  mutate(p.adj = -log10(p.adj), p.adj = ifelse(type == "Down", p.adj * -1, p.adj)) %>%
  ggplot(aes(x = term_name, y = Treatment)) +
  geom_point(aes(size = intersection_size, fill = p.adj), shape = 21) +
  xlab("") + ylab("") + ggtitle("Top Enriched Gene Sets") +
  lgd +
  theme(axis.text = element_text(size = 10), 
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1, size = 8.5), 
        panel.background = element_rect(linewidth = 1, fill = "transparent")) +
  scale_size_continuous(name = "Number of \nOverlapping Genes", 
                        breaks = seq(0, max(gseaTab$intersection_size), by = 50), 
                        limits = c(1, 200), range = c(0.5,7.5)) +
  scale_fill_gradient2(low = "#1976D2", high = "#F44336", limits=c(-4, 4), oob=squish, 
                       name = "-Log10(p.adj) \nwith Direction")
  
p

ggsave(plot = p, filename = paste0(opt$plot, "top_genesets_dotplot.png"), 
       height = 19.5, width = 50, units = "cm")

Output session info

## R version 4.3.2 (2023-10-31)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
##  [1] parallel  grid      stats4    stats     graphics  grDevices datasets 
##  [8] utils     methods   base     
## 
## other attached packages:
##  [1] sna_2.7-2                   statnet.common_4.9.0       
##  [3] scales_1.3.0                numbat_1.4.0               
##  [5] Matrix_1.6-1.1              airr_1.5.0                 
##  [7] ggsci_3.2.0                 BiocParallel_1.36.0        
##  [9] ggnet_0.1.0                 network_1.18.2             
## [11] tftargets_1.3               ggrepel_0.9.5              
## [13] readxl_1.4.3                edgeR_4.0.16               
## [15] limma_3.58.1                circlize_0.4.16            
## [17] ComplexHeatmap_2.18.0       gridExtra_2.3              
## [19] gprofiler2_0.2.3            scater_1.30.1              
## [21] scran_1.30.2                scuttle_1.12.0             
## [23] SingleCellExperiment_1.24.0 SummarizedExperiment_1.32.0
## [25] Biobase_2.62.0              GenomicRanges_1.54.1       
## [27] GenomeInfoDb_1.38.8         IRanges_2.36.0             
## [29] S4Vectors_0.40.2            BiocGenerics_0.48.1        
## [31] MatrixGenerics_1.14.0       matrixStats_1.3.0          
## [33] lubridate_1.9.3             forcats_1.0.0              
## [35] stringr_1.5.1               dplyr_1.1.4                
## [37] purrr_1.0.2                 readr_2.1.5                
## [39] tidyr_1.3.1                 tibble_3.2.1               
## [41] ggplot2_3.5.1               tidyverse_2.0.0            
## 
## loaded via a namespace (and not attached):
##   [1] bitops_1.0-7              ggplotify_0.1.2          
##   [3] cellranger_1.1.0          polyclip_1.10-7          
##   [5] lifecycle_1.0.4           doParallel_1.0.17        
##   [7] lattice_0.21-9            MASS_7.3-60              
##   [9] magrittr_2.0.3            plotly_4.10.4            
##  [11] sass_0.4.9                rmarkdown_2.27           
##  [13] jquerylib_0.1.4           yaml_2.3.9               
##  [15] metapod_1.10.1            RColorBrewer_1.1-3       
##  [17] abind_1.4-5               zlibbioc_1.48.2          
##  [19] quadprog_1.5-8            hahmmr_1.0.0             
##  [21] ggraph_2.2.1              RCurl_1.98-1.14          
##  [23] yulab.utils_0.1.5         tweenr_2.0.3             
##  [25] GenomeInfoDbData_1.2.11   irlba_2.3.5.1            
##  [27] tidytree_0.4.6            dqrng_0.4.1              
##  [29] DelayedMatrixStats_1.24.0 codetools_0.2-19         
##  [31] DelayedArray_0.28.0       ggforce_0.4.2            
##  [33] tidyselect_1.2.1          shape_1.4.6.1            
##  [35] aplot_0.2.3               farver_2.1.2             
##  [37] ScaledMatrix_1.10.0       viridis_0.6.5            
##  [39] jsonlite_1.8.8            GetoptLong_1.0.5         
##  [41] BiocNeighbors_1.20.2      tidygraph_1.3.1          
##  [43] iterators_1.0.14          systemfonts_1.1.0        
##  [45] foreach_1.5.2             tools_4.3.2              
##  [47] ragg_1.3.2                treeio_1.29.0.002        
##  [49] Rcpp_1.0.12               glue_1.7.0               
##  [51] SparseArray_1.2.4         xfun_0.45                
##  [53] withr_3.0.0               fastmap_1.2.0            
##  [55] bluster_1.12.0            fansi_1.0.6              
##  [57] digest_0.6.36             rsvd_1.0.5               
##  [59] parallelDist_0.2.6        timechange_0.3.0         
##  [61] R6_2.5.1                  gridGraphics_0.5-1       
##  [63] textshaping_0.4.0         colorspace_2.1-0         
##  [65] RhpcBLASctl_0.23-42       utf8_1.2.4               
##  [67] generics_0.1.3            renv_1.0.7               
##  [69] data.table_1.15.4         graphlayouts_1.1.1       
##  [71] httr_1.4.7                htmlwidgets_1.6.4        
##  [73] S4Arrays_1.2.1            pkgconfig_2.0.3          
##  [75] gtable_0.3.5              XVector_0.42.0           
##  [77] htmltools_0.5.8.1         clue_0.3-65              
##  [79] png_0.1-8                 ggfun_0.1.5              
##  [81] knitr_1.48                rstudioapi_0.16.0        
##  [83] tzdb_0.4.0                rjson_0.2.21             
##  [85] coda_0.19-4.1             nlme_3.1-163             
##  [87] cachem_1.1.0              GlobalOptions_0.1.2      
##  [89] vipor_0.4.7               pillar_1.9.0             
##  [91] logger_0.3.0              vctrs_0.6.5              
##  [93] BiocSingular_1.18.0       beachmat_2.18.1          
##  [95] cluster_2.1.4             beeswarm_0.4.0           
##  [97] evaluate_0.24.0           cli_3.6.3                
##  [99] locfit_1.5-9.10           compiler_4.3.2           
## [101] rlang_1.1.4               crayon_1.5.3             
## [103] labeling_0.4.3            fs_1.6.4                 
## [105] ggbeeswarm_0.7.2          stringi_1.8.4            
## [107] viridisLite_0.4.2         munsell_0.5.1            
## [109] lazyeval_0.2.2            hms_1.1.3                
## [111] patchwork_1.2.0           sparseMatrixStats_1.14.0 
## [113] statmod_1.5.0             highr_0.11               
## [115] igraph_2.0.3              memoise_2.0.1            
## [117] scistreer_1.2.0           RcppParallel_5.1.8       
## [119] bslib_0.7.0               phangorn_2.11.1          
## [121] ggtree_3.10.1             fastmatch_1.1-4          
## [123] ape_5.8